Books

• Jean Claude, Morphometry with R
  • Online through ETHZ
  • Buy it
• John C. Russ, “The Image Processing Handbook”,(Boca Raton, CRC Press)
  • Available online within domain ethz.ch (or proxy.ethz.ch / public VPN)
• J. Weickert, Visualization and Processing of Tensor Fields
  • Online

Papers / Sites

• Voronoi Tesselations

  • Ghosh, S. (1997). Tessellation-based computational methods for the characterization and analysis of heterogeneous microstructures. Composites Science and Technology, 57(9-10), 1187–1210
  • Wolfram Explanation

• Self-Avoiding / Nearest Neighbor

  • Schwarz, H., & Exner, H. E. (1983). The characterization of the arrangement of feature centroids in planes and volumes. Journal of Microscopy, 129(2), 155–169.
  • Kubitscheck, U. et al. (1996). Single nuclear pores visualized by confocal microscopy and image processing. Biophysical Journal, 70(5), 2067–77.

• Alignment / Distribution Tensor

  • Mader, K. et al (2013). A quantitative framework for the 3D characterization of the osteocyte lacunar system. Bone, 57(1), 142–154

• Two point correlation

  • Dinis, L., et. al. (2007). Analysis of 3D solids using the natural neighbour radial point interpolation method. Computer Methods in Applied Mechanics and Engineering, 196(13-16)

• Motivation (Why and How?)
• Local Environment
  • Neighbors
  • Voronoi Tesselation
  • Distribution Tensor

• Global Enviroment
  • Alignment
  • Self-Avoidance
  • Two Point Correlation Function

With all of these images, the first step is always to understand exactly what we are trying to learn from our images.

1. We want to know how many cells are alive

   • Maybe small cells are dead and larger cells are alive \( \rightarrow \) examine the volume distribution
   • Maybe living cells are round and dead cells are really spiky and pointy \( \rightarrow \) examine anisotropy

2. We want to know where the cells are alive or most densely packed

   • We can visually inspect the sample (maybe even color by volume)
   • We can examine the raw positions (x,y,z) but what does that really tell us?
   • We can make boxes and count the cells inside each one
   • How do we compare two regions in the same sample or even two samples?

1. We want to know how the cells are communicating

   • Maybe physically connected cells (touching) are communicating \( \rightarrow \) watershed
   • Maybe cells oriented the same direction are communicating \( \rightarrow \) average? orientation
   • Maybe cells which are close enough are communicating \( \rightarrow \) ?
   • Maybe cells form hub and spoke networks \( \rightarrow \) ?

2. We want to know how the cells are nourished

   • Maybe closely packed cells are better nourished \( \rightarrow \) count cells in a box?
   • Maybe cells are oriented around canals which supply them \( \rightarrow \) ?

1. A way for counting cells in a region and estimating density without creating arbitrary boxes
2. A way for finding out how many cells are near a given cell, it's nearest neighbors
3. A way for quantifying how far apart cells are and then comparing different regions within a sample
4. A way for quantifying and comparing orientations
   • the mean of 0\( ^\circ \) and 180\( ^\circ \) = 90\( ^\circ \)
   • the distance between -180\( ^\circ \) and 179\( ^\circ \) is 359\( ^\circ \)
   • since we have not defined a tip or head, 0\( ^\circ \) and 180\( ^\circ \) are actually the same

What would be really great?

A tool which could be adapted to answering a large variety of problems

• multiple types of structures
• multiple phases

With most imaging techniques and sample types, the task of measurement itself impacts the sample.

• Even techniques like X-ray tomography which claim to be non-destructive still impart significant to lethal doses of X-ray radition for high resolution imaging
• Electron microscopy, auto-tome-based methods, histology are all markedly more destructive and make longitudinal studies impossible
• Even when such measurements are possible
  • Registration can be a difficult task and introduce artifacts

Why is this important?

• techniques which allow us to compare different samples of the same type.
• are sensitive to common transformations
  • Sample B after the treatment looks like Sample A stretched to be 2x larger
  • The volume fraction at the center is higher than the edges but organization remains the same

\[ \downarrow \]

\[ \cdots \]

So if we want to know the the mean or standard deviations of the position or orientations we can analyze them easily.

• But what if we want more or other information?

One of the first metrics to examine with distribution is density \( \rightarrow \) how many objects in a given region or volume.

It is deceptively easy to calculate involving the ratio of the number of objects divided by the volume.

It doesn't tell us much, many very different systems with the same density and what if we want the density of a single point? Does that even make sense?

Definition

Oxford American \( \rightarrow \) be situated next to or very near to (another)

• Does not sound very scientific
• How close?
  • Touching, closer than anything else?

Nearest Neighbor (distance)

Given a set of objects with centroids at \[ \textbf{P}=\begin{bmatrix} \vec{x}_0,\vec{x}_1,\cdots,\vec{x}_i \end{bmatrix} \]

We can define the nearest neighbor as the position of the object in our set which is closest

\[ \vec{\textrm{NN}}(\vec{y}) = \textrm{argmin}(||\vec{y}-\vec{x}|| \forall \vec{x} \in \textbf{P}-\vec{y}) \]

We define the distance as the Euclidean distance from the current point to that point, and the angle as the

\[ \textrm{NND}(\vec{y}) = \textrm{min}(||\vec{y}-\vec{x}|| \forall \vec{x} \in \textbf{P}-\vec{y}) \] \[ \textrm{NN}\theta(\vec{y}) = \tan^{-1}\frac{(\vec{\textrm{NN}}-\vec{y})\cdot \vec{j}}{(\vec{\textrm{NN}}-\vec{y})\cdot \vec{i}} \]

So examining a simple starting system like a grid, we already start running into issues.

• In a perfect grid like structure each object has 4 equidistant neighbors (6 in 3D)
• Which one is closest?

We thus add an additional clause (only relevant for simulated data) where if there are multiple equidistant neighbors, a nearest is chosen randomly

This ensures when we examine the orientation distribution (NN\( \theta \)) of the neighbors it is evenly distributed

For the rest of these sections we will repeatedly use several simple in-silico systems to test our methods and try to better understand the kind of results we obtain from them.

• Compression

  • The most simple system simply involves a scaling in every direction by \( \alpha \)
  • \( \alpha<1 \) the system is compressed
  • \( \alpha>1 \) the system is expanded

  \[ \begin{bmatrix} x^\prime \\ y^\prime \end{bmatrix} = \alpha \begin{bmatrix} x \\ y \end{bmatrix} \]

• Stretch
  • A non-evenly distributed system with a parameter \( \alpha \) controlling if objects bunch near the edges or the center. A maximum distance \( m \) is defined as the magnitude of the largest offset in x or y
  • + Same total volume just arranged differently

\[ \begin{bmatrix} x^\prime \\ y^\prime \end{bmatrix} = \begin{bmatrix} \textrm{sign}(x) \left(\frac{|x|}{m}\right)^\alpha m \\ \textrm{sign}(y) \left(\frac{|y|}{m}\right)^\alpha m \end{bmatrix} \]

• Swirl
  • A transformation where the points are rotated more based on how far away they are from the center and the slope of the swirl (\( \alpha \)),

\[ \theta (x,y) = \alpha \sqrt{x^2+y^2} \]

\[ \begin{bmatrix} x^\prime \\ y^\prime \end{bmatrix} = \begin{bmatrix} \cos\theta(x,y) & -\sin\theta(x,y) \\ \sin\theta(x,y) & \cos\theta(x,y) \end{bmatrix} \begin{bmatrix} x \\ y \end{bmatrix} \]

Using a uniform grid of points as a starting point has a strong influence on the results. A better approach is to use a randomly distributed series of points

• resembles real data much better
• avoids these symmetry problems
  • \( \epsilon \) sized edges or overlaps
  • identical distances to nearby objects

Voronoi tesselation is a method for partitioning a space based on points. The basic idea is that each point \( \vec{p} \) is assigned a region \( \textbf{R} \) consisting of points which are closer to \( \vec{p} \) than any of the other points. Below the diagram is shown in a dashed line for the points shown as small circles.

We call the area of a region (\( \textbf{R} \)) around point \( \vec{p} \) its territory.

The grid on the random system, shows much more diversity in territory area.

A parallel or dual idea where triangles are used and each triangle is created such that the circle which encloses it contains no other points. The triangulation makes the neighbors explicit since connected points in the triangulation correspond to points in our tesselation which share an edge (or face in 3D)

We define the number of connections each point \( \vec{p} \) has the Neighbor Count or Delaunay Neighbor Count.

The triangulation on a random system has a much higher diversity in neighbor count

Compression

Tension

Low Shear

High Shear

Low Stretch

Highly Stretched System

Low Swirl System

High Swirl System

Compression

Stretch

Shear

Swirl

We have introduced a number of “operations” we can perform on our objects to change their positions

• compression
• stretching
• shearing
• swirling

We have introduced a number of metrics to characterize our images

• Nearest Neighbor distance
• Nearest Neighbor angle
• Delaunay Neigbhor count
• Territory Area (Volume)

A single random systems is useful

• but in order to have a reasonable understanding of the behavior of a system we need to sample many of them.

Understand metrics as a random system + a known transformation

• We take mean values
• Also coefficient of variation (CV) values since they are “scale-free”

In imaging science we always end up with lots of data, the tricky part is understanding the results that come out. With this simulation-based approach

• we generate completely random data
• apply a known transformation to it \( \mathcal{F} \)
• quantify the results

We can then take this knowledge and use it to interpret observed data as transformations on an initially random system. We try and find the rules used to produce the sample

Examples

1. Cell distribution in bone

   • Cell position appears random
   • Metrics \( \neq \) Random Statistics
   • Cells are consistently self-avoiding or stretched

2. Egg-shell Pores

   • Pores in rock / egg shell appear random
   • Metrics \( \neq \) Random Statistics
   • Pores are also self-avoiding

From the nearest neighbor distance metric, we can create a scale-free version of the metric which we call self-avoiding coefficient or grouping.

The metric is the ratio of

• observed nearest neighbor distance \( \textrm{NND} \)
• the expected mean nearest neighbor distance (\( r_0 \)) for a random point distribution (Poisson Point Process) with the same number of points (N.Obj) per volume (Total.Volume). \[ r_0=\sqrt[3]{\frac{\textrm{Total.Volume}}{2\pi \textrm{ N.Obj}}} \]

Using the territory we defined earlier (Region area/volume) we can simplify the definition to

\[ r_0=\sqrt[3]{\frac{\bar{\textrm{Ter}}}{2\pi}} \]

\[ \textrm{SAC} = \frac{\textrm{NND}}{\sqrt[3]{\frac{\bar{\textrm{Ter}}}{2\pi}}} \]

So the information we have is 3D why are we taking single metrics (distance, angle, volume) to quantify it.

• Shouldn't we use 3D metrics with 3D data?
• Just like the shape tensor we covered before, we can define a distribution tensor to characterize the shape of the distribution.
• The major difference instead of constituting voxels we use edges
  • an edge is defined from the Delaunay triangluation
  • it connects two neighboring bubbles together

• We can calculate distribution for a single bubble by taking all edges that touch the object or from a region by finding all edges inside that region

We start off by calculating the covariance matrix from the list of edges \( \vec{v}_{ij} \) in a given volume \( \mathcal{V} \)

\[ \vec{v}_{ij} = \vec{\textrm{COV}}(i)-\vec{\textrm{COV}}(j) \]

\[ COV(\mathcal{V}) = \frac{1}{N} \sum_{\forall \textrm{COM}(i) \in \mathcal{V}} \begin{bmatrix} \vec{v}_x\vec{v}_x & \vec{v}_x\vec{v}_y & \vec{v}_x\vec{v}_z\\ \vec{v}_y\vec{v}_x & \vec{v}_y\vec{v}_y & \vec{v}_y\vec{v}_z\\ \vec{v}_z\vec{v}_x & \vec{v}_z\vec{v}_y & \vec{v}_z\vec{v}_z \end{bmatrix} \]

Visual example

• Tensor represents the average spacing between objects in each direction \( \approx \) thickness of background.
• Its interpretation is more difficult since it doesn't represent a real object

From this tensor we can define an anisotropy in the same manner as we defined for shapes. The anisotropy defined as before \[ Aiso = \frac{\text{Longest Side}-\text{Shortest Side}}{\text{Longest Side}} \]

• An isotropic distribution tensor indicates the spacing between objects is approximately the same in every direction
  • Does not mean a grid, organized, or evenly spaced!
  • Just the same in every direction
• Anisotropic distribution tensor means the spacing is smaller in one direction than the others
  • Can be regular or grid-like
  • Just closer on one direction

From this tensor we can also define oblateness in the same manner as we defined for shapes. The oblateness is also defined as before as a type of anisotropy

\[ \textrm{Ob} = 2\frac{\lambda_{2}-\lambda_{1}}{\lambda_{3}-\lambda_{1}}-1 \]

• Indicates the pancakeness of the distribution
• Oblate means it is very closely spaced in 2 directions and far in the other
  • strand-like structure
• Prolate means it is close in one direction and far in the other 2
  • sheet-like structure

The shape tensor provides for each object 3 possible orientations (each of the eigenvectors).

As we saw in earlier plots it is difficult to express this information in point or graph format

We can again take advantage of the versatility of a tensor representation for our data and use an alignment tensor.

• The same as all of the tensors we have introduced
• Except we use the primary orientation instead of voxel positions (shape) or edges (distribution)
• Similar to distribution it is calculated for a volume (\( \mathcal{V} \)) or region and is meaningless for a single object

\[ \vec{v}_{i} = \vec{\Lambda_1}(i) \]

\[ COV = \frac{1}{N} \sum_{\forall \textrm{COM}(i) \in \mathcal{V}} \begin{bmatrix} \vec{v}_x\vec{v}_x & \vec{v}_x\vec{v}_y & \vec{v}_x\vec{v}_z\\ \vec{v}_y\vec{v}_x & \vec{v}_y\vec{v}_y & \vec{v}_y\vec{v}_z\\ \vec{v}_z\vec{v}_x & \vec{v}_z\vec{v}_y & \vec{v}_z\vec{v}_z \end{bmatrix} \]

K-Means

K-Means can also be used to classify the point-space itself after shape analysis. It is even better suited than for images because while most images are only 2 or 3D the shape vector space can be 50-60 dimensional and is inherently much more difficult to visualize.

2 Point Correlation Functions

For a wider class of analysis of spatial distribution, there exist a class a functions called N-point Correlation Functions

• given a point of type \( A \) at \( \vec{x} \)

• what is the probability of \( B \) at \( \vec{x}+\vec{r} \)

• How is it useful

  • A simple analysis can be used to see the spacing of objects (peaks in the \( F(\vec{r}) \) function, \( A \) and \( B \) are the same phase).
  • Where are cells located in reference to canals (peaks in \( F(\vec{r}) \) function, \( A \) is a point inside a canal \( B \) is the cells)